Spin coupling around a carbon atom vacancy in graphene 
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We investigate the details of the electronic structure in the neighborhoods of a carbon atom 
vacancy in graphene by employing magnetization-constrained density-functional theory on periodic 
slabs, and spin-exact, multi-reference, second-order perturbation theory on a finite cluster. The 
picture that emerges is that of two local magnetic moments (one 7r-like and one cr-like) decoupled 
from the ty band and coupled to each other. We find that the ground state is a triplet with a planar 
equilibrium geometry where an apical C atom opposes a pentagonal ring. This state lies ^0.2 eV 
lower in energy than the open-shell singlet with one spin nipped, which is a bistable system with 
two equivalent equilibrium lattice configurations (for the apical C atom above or below the lattice 
plane) and a barrier ^0.1 eV high separating them. Accordingly, a bare carbon-atom vacancy is 
predicted to be a spin-one paramagnetic species, but spin-half paramagnetism can be accommodated 
if binding to foreign species, ripples, coupling to a substrate, or doping are taken into account. 



I. INTRODUCTION 

Magnetism in graphene is a fascinating and highly 
controversial mattepR Early reports on ferromagnetic 
ordering in graphite and graphene^^ have been ques- 
tioned in the light of the ubiquitous presence of mag- 
netic contaminants, and measurements under care- 
fully controlled conditions showed that graphene, like 
graphite, is strongly diamagnetic with a weak para- 
magnetic contribution from adatoms and/or carbon 
atom vacancies 6 . Simple adsorbates such as fluo- 
rine and missing carbon atoms have been shown to 
provide a spin- 1/2 paramagnetic response^, though 
spin-1 paramagnetism has been reported upon N + 
irradiation 8 . Direct evidences of the presence of 
adatom- and vacancy-related magnetic moments have 
been observed in pure spin transport measurements 
on graphene spin valves 9 . Likewise, signatures of 
Kondo effect have been observed in charge-transport 
measurements at low temperatures on irradiated 
graphene, and high Kondo temperatures reported for 
finite and zero electron density^, at odds with the 
above magnetometry measurements 7 . 

From a theoretical perspective, perfect bipartite 
systems support a number of zero-energy "midgap" 
states which is greater or equal than the sublattice 
imbalance \ua — tib\, where n^,n g are the number 
of sites in the two sublattice d 11 * 12 *. When imbalance 
results from isolated missing p z orbitals (e.g. for low 
concentrations of covalently bound adatoms or vacan- 
cies) these states decay slowly (~ 1/r) from t he de- 
fects and localize on the locally majority sites 13 14 , 
as confirmed by experiments^. Thus, these defects 
form quasi- localized tt moments, which couple either 
ferromagnetically or antiferromagnetically depending 
on their lattice position. In fact, with local interac- 
tions only, at charge neutrality (half-filling) the spin 
state of the system exactly matches the sublattice 
imbalance^, S = — ng|/2, and thus coupling is 
ferromagnetic for defects in the same sublattice and 
antiferromagnetic otherwise, a result which also fol- 



lows by analyzing RRKY interactions in graphene^. 
Within the same assumptions (perfect electron-hole 
symmetry, local interactions only) coupling between 
7r moments and conduction states has been inves- 
tigated beyond mean-field approaches and found to 
be ferromagnetic^, thereby confirming that simple 
adatoms covalently bound to the substrate (e.g. H, 
F species) behave as spin-1/2 localized moments. In 
turn, this also affects chemical prop erties and favors 
formation of dimers of balanced typ d 19 * 20 -! 

This simple picture has to be revised for a car- 
bon atom vacancy where, in addition to the above 
7r midgap state, three a orbitals are left singly oc- 
cupied upon vacancy formation, and a structural in- 
stability (Jahn- Teller distortion) arises which breaks 
electron-hole symmetry, even if nearest neighbors in- 
teractions only are retained. The ensuing lattice 
re- arrangement leaves two unpaired electrons, and 
a magnetic moment in the range 2.0 — 1.0 \±b has 
been found by (ensemble) density functional theory 
(DFT) calculations 2 ^ 2 ^, with a clear tendency to 
1.0 He (and vanishing dependence of the energy on 
the magnetization) in the low-density limit which sig- 
nals the absence of any order at experimentally rel- 
evant concentrations 2 ^. Apart from the possible role 
of electron-hole symmetry, this apparently contrasts 
with the situation described above for a tt moment 
only, since the presence of an additional unpaired 
electron in a very localized a orbital is expected to 
give either a triplet or a singlet state. Recent experi- 
ments have shown that the spin— 1/2 paramagnetism 
of missing carbon atoms has two contributions 2 ^, from 
a and tt states respectively, and one of them can 
be quenched upon molecular doping and possibly by 
means of the electric field effect! 2 - 7 -'. This result re- 
quires that, at least under the experimental conditions 
of Ref. [27l the unpaired electrons around a vacancy 
negligibly interact with each other, in contrast with 
early reports on spin-1 paramagnetism of irradiated 
graphene samples 8 ^ Proper consideration of a states, 
and their possible hybridization with tt states when 
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the substrate is no longer locally planar, e.g. because 
of ripples or interaction with a substrate 2 ^, has led 
to reconsidering the issue of the interaction between 
the localized magnetic moments and the conduction 
electron d 29 * 30 *, though the most recent experiments^ 
seem to rule out this possibility. 

Here, in order to help shed light on the above issues 
we consider in detail the electronic structure around 
a carbon atom vacancy in graphene. We employ both 
conventional DFT methods in periodic models and 
accurate, spin-exact quantum chemistry methods in 
a finite cluster, to investigate the coupling of the two 
electrons left unpaired after vacancy formation and re- 
construction. We analyze several substrate geometries 
close to the equilibrium one, and focus in particular on 
the out-of-plane movement of the carbon atom where 
most of the unpaired electron density resides. Our 
results show that the triplet is the ground-state and 
has a planar equilibrium geometry, while the singlet 
-which lies ~ 0.2 eV above it- gets easily stabilized by 
an out-of-plane lattice distortion. Hence, both spin-1 
and spin- 1/2 paramegnetism may in principle arise in 
irradiated graphene, depending on local interactions, 
curvature, etc. of the graphene sheet, in addition to 
doping or chemical interactions with foreign species. 

The manuscript is organized as follows: Section [IT] 
outlines the important Jahn- Teller distortion occur- 
ring in the system and Section [ITT] gives the details of 
the electronic structure methods adopted in this work. 
Results are given in Section [IV| and discussed in Sec- 
tion V, and Section VI summarizes and concludes. 



II. JAHN-TELLER DISTORTION 

As mentioned above, formation of a carbon atom 
vacancy gives rise to localized states around the va- 
cancy, namely one tt (semilocalized) "midgap" state 
and three dangling orbitals in the a network which 
result from breaking the sp 2 bonds which hold the 
carbon atom in place. In the local Dsh point sym- 
metry group which is appropriate to discuss proper 
and pseudo Jahn- Teller distortions, the first belongs 
to a 2 symmetry species, and the latter span a[ + e' 
irreducible representations, a[ being lowest in en- 
ergy since it contains a purely bonding combina- 
tion of a orbitals. Hence, the lowest energy "scenar- 
ios" for the many-body electronic state can be ob- 
tained by distributing two electrons in the e' and a 2 f 
states, i.e. starting from configurations of the type 

..(a / 1 ) 2 (e / ) ni ( a 2 / ) n2 with n i + n< 2 = 2 - Among these, 
the one with m = 712 = 1 is expected to be low- 
est in energy and gives rise to many-body states of 
E" symmetry for both the parallel and antiparallel 
alignment [The remaining possibilities with two elec- 
trons in the same set of states are pushed up in en- 
ergy by a larger Coloumb repulsion and have sym- 
metries 1 A' 1 + 3 A! 2 + 1 E' for m = 2 and 1 A[ for 




Figure 1: Optimized structure of a carbon vacancy in a 
6x6 unit cell, with the apical carbon marked by an arrow. 



ri2 = 2]. Thus, the ground-state is doubly degener- 
ate for both spin alignments and undergoes (proper 
or pseudo) Jahn- Teller distortion. The latter oc- 
curs because of coupling with in-plane e' vibrations 
(\E"f = [E'} 2 = A' + E') which distort the symmetric 
arrangement of the carbon atoms around the vacancy. 
This is a standard E<&e problem which is described by 
the so-called tricor n whe n such vibrations are included 
up to second-orde j 31 * 32 l 



Accordingly, a distorted ge- 
ometry with a "pentagonal" ring and a "apical" carbon 
atom opposite to it (see Fig{T]), can be predicted to be 
the equilibrium configuration (three-fold degenerate), 
as indeed found in several previous investigations^"^. 

Out of plane, E" vibrations do not lift degeneracy at 
first-order, but may affect energetics at higher orders, 
especially if coupling to the low-lying excited states 
is considered^. This is particularly important here 
since such distortions are qualitatively different for the 
states depending on whether the spins are par- 

allel or antiparallel to each other. This is evident for 
the out-of-plane movement of the apical carbon atom 
in the distorted configuration: the a and the 7r states 
may hybridize to some extent (and gain energy from 
double filling) if the two electrons couple at low-spin, 
otherwise they require extra energy to reduce their 
overlap. As a consequence, the planar structure is ex- 
pected to be stable in the triplet state only; in the sin- 
glet, a non-planar configuration with the apical carbon 
atom slightly above (or below) the surface plane ap- 
pears to be more stable. Because of that, the relative 
stability of the two spin states is geometry-dependent 
and its analysis requires at least investigating the out- 
of-plane movement of the "apical" carbon atom. This 



is described in Section IIVI after Section III has intro- 



duced the electronic structures methods and setups 
adopted, along with the structural models chosen to 
investigate the vacancy. 



III. METHODS AND MODELS 

Electronic structure calculations were performed 
at different correlation levels for different structural 
models. Periodic arrangements of vacancies in large 
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Table I: Results of full structural relaxation without con- 
straints on the magnetization, for a vacancy in several 
n x n supercells. AE is the energy separation between 
the metastable non-planar (C s ) structure and the planar 
(C2v) minimum, M np is the magnetization of the former 
and M p that of the latter. Also reported the length of 
newly formed CC bond closing the pentagon (d p cc and 
dc C f° r planar and non-planar geometries, respectively) 
and the height he of the apical carbon atom in the non- 
planar configuration. 

unit cells were investigated with standard, plane- wave 
based density functional theory calculations, whereas 
a finite-size (cluster) model was judiciously selected 
and studied with correlated wavefunction methods de- 
scribed below. 




Figure 2: The molecular model adopted for the wavefunc- 
tion calculations, along with isosurfaces of singly occupied 
a (left) and tt -midgap (right) orbitals for \<f>\ — 0.015 
A -3 / 2 , as obtained at the restricted open-shell Hartree- 
Fock level (S = 1) for the minimum structure. 



for the minimum are in good agreement with those 
found in previous studies^^l. 

On the basis of these results we concluded that a 
6x6 super cell with a 6 x 6 x 1 k— point mesh was 
a good compromise between the need of reducing in- 
teraction between periodic images and computational 
manageability. Therefore, further investigations were 
performed with this setup. 



B. A finite-size model 



A. Periodic models 

Periodic models were studied with plane- wave DFT 
as impl emen ted in the Vienna ab initio package suite 
(VASP; 33 l 34 l The exchange-correlation effects were 
treated with the Perdew-Burke-Ernzerhof (PBE) 3 -^ 3 -^ 
functional within the generalized gradient approxima- 
tion (GGA), in the spin-polarized framework. Kohn- 
Sham orbitals were expended on a plane- wave set lim- 
ited to a 500 eV energy cutoff and core electrons were 
frozen and replaced by projector- augmented wave 
(PAW) pseudopotentials^ 3 -^. Several n x n graphene 
supercells with a 20 A vacuum were initially consid- 
ered to model the defective system, from n = 2 to 
n = 10, by using T centered k— point meshes rang- 
ing from 15 x 15 x 1 (for n = 2) to 3 x 3 x 1 for 
n = 6 — 10. The structure of a vacancy in such cells 
was fully optimized without constraints on the magne- 
tization and gave a Jahn- Teller distorted planar min- 
imum, with a C — C bond length in the pentagon 
of about 2 A, i.e. smaller than the graphene lattice 
constant a — 2.46 A but much larger than a typical 
(single) C — C bond (1.54 A). Additionally, starting 
with a low-magnetization guess, we invariably found 
for n > 4 a metastable spin-polarized solution with 
vanishing magnetization and a non-planar configura- 
tion. Its energy separation AE from the planar min- 
imum in the same supercell, along with the resulting 
magnetizations M and the main geometrical parame- 
ters of the two structures are given in Table [TJ Results 



The semi-localized character of the electronic states 
induced by the vacancy makes it realistic the study 
of finite-size models where accurate, correlated wave- 
function calculations are feasible. The size (and 
shape) of the cluster has to be large enough to min- 
imize the effects that edges have on the details of 
the electronic structure, and small enough that com- 
plex many-body wavefunctions are yet tractable. To 
this end we considered a reasonably sized Polycyclic 
Aromatic Hydrocarbon (PAH) a carbon 

cluster with a central vacant site which is hydrogen- 
terminated at the edges (Fig. Its actual shape was 
chosen -following the line of reasoning of Ref. l39l with 
the help of Tight-Binding (TB) calculations in such a 
way to limit the edge localization which does inter- 
fere with the defect-induced states at the Fermi level. 
In the chosen structure, edge states were found suffi- 
ciently far in energy from the vacancy-induced states 
(both at the TB and at the Hartree-Fock (HF) level) 
to make us confident that the resulting energetics ac- 
curately describes the vacancy. Cluster geometries 
were selected with a 'cut-out process' starting from the 
above mentioned 6x6 supercells, and adding hydrogen 
atoms to the undercoordinated edge C atoms, without 
further geometry refinement. In this way, comparison 
between the periodic and the cluster model with the 
same local arrangement close to the vacancy was pos- 
sible. 

Accurate results on the finite model were achieved 
through all-electron, correlated wavefunction calcula- 
tions based on atom-centered basis-sets of the correla- 
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tion consistent type (cc-pVDZ). Energy was obtained 
with the help of the MOLPRO suite of codes 40 by cor- 
recting to second order in perturbation theory a 'refer- 
ence' wavefunction of the Complete Active Space Self- 
Consistent Field (CASSCF) type, according to what 
is known as CASPT^EU. The CASSCF(n,ra) wave- 
function is a multi-determinant wavefunction contain- 
ing all possible excitations of n 'active' electrons in 
m 'active' orbitals, where all the orbit als a nd e xpan- 
sion coefficients are variationally optimized 43 * 44 *. For 
our purposes, we started with a minimal active space 
containing the a and tt orbitals localized on the va- 
cancy and the two electrons occupying them at the 
HF level, and enlarged it by including two further tt 
orbitals (one below and one above the Fermi level), i.e. 
CAS (4,4). We fully optimized the active orbitals and 
the thirty doubly-occupied orbitals higher in energy, 
and kept the lowest-lying (doubly occupied) orbitals 
frozen at the Hartree-Fock level. With the optimized 
CASSCF wavefunctions at hand, dynamic correlation 
was introduced by including perturbatively the effect 
of single and double excitations out of the configura- 
tions contained in the selected CAS 'reference' space. 



IV. RESULTS 

A. Periodic calculations 

As shown in Tabjl] all the optimized structures 
show an appreciable magnetization, in agreement with 
previous results obtained for similar systemP^H, 
and this is significantly different for the two struc- 
tures. The two geometries thus refer to the differ- 
ent electronic states -the high and the low spin state, 
respectively- which originate from the different ways 
the unpaired electrons localized on the vacancy align. 
Unfortunately, even assuming that the KS wavefunc- 
tion correctly describes the spin state of the sys- 
tem, proper assignment is not possible with DFT. 
This is because (i) allowance of fractional occupation 
of single-particle Kohn-Sham states, as occurs with 
ensemble-DFT, typically prevents the appearance of 
pure states even when they are in principle possible 
at T — K, and (ii) pure but spin-polarized (unre- 
stricted) solutions do not have a definite spin value 

In order to get rid of these problems, and to mimic 
as much as possible the expected electron configu- 
rations, we relied on magnetization-constrained DFT 
calculations on the 6x6 supercell, setting the (pro- 
jection of the) magnetic moment to two (zero) Bohr 
magnetons for the triplet (singlet) case. Full struc- 
tural optimizations were then performed for different 
out-of-plane displacements h of the apical carbon, for 
each 'spin' state, to investigate how these states evolve 
out of the plane. 

The results of such calculations are shown in Fig. 




Figure 3: Magnetization-constrained energies as functions 
of the height h of the apical carbon atom above the surface. 
Filled and empty symbols for M = 2,0 /is, referenced 
to the minimum of the M — 2\ib case. Also reported 
as a thick line the results of magnetization-unconstrained 
calculations, referenced to their minimum, and the corre- 
sponding magnetization (dashed line, right scale). 



[3J referenced to the planar configuration in the triplet 
state, along with the magnetization-unconstrained 
curve referenced to its minimum (which is only 29 
meV below that of the constrained triplet curve). It is 
clear that the latter is a mixture of the two electronic 
states, with the triplet prevailing for h ~ and the 
singlet dominant for h ^> 0. 

From Fig. |3j we also see that the global minimum 
belongs to the triplet curve and has a flat geometry 
(of symmetry C^). The singlet curve instead shows 
two equivalent minima (of symmetry C s ) for the car- 
bon atom above and below the plane, respectively, 
~ 0.4 A away from the surface. The latter thus repre- 
sents a bi-stable system that crosses the triplet when 
the carbon atom moves out by about ~ 0.5 A, but 
is otherwise higher in energy. The energy difference 
between the singlet and triplet minima is ~ 0.18 eV, 
thus significantly larger than the (unconstrained) AE 
for the same 6x6 supercell reported in Table [TJ which 
referred to 'mixed' electronic states. The singlet min- 
ima are separated by a barrier ~ 0.2 eV high which 
lies ~ 0.4 eV above the triplet. This estimate will be 
refined in the next section on the basis of more accu- 
rate electronic structure calculations. Before leaving 
this section we only stress that the curves reported 
in Fig. [3] refer to a full structural relaxation (in the 
given electronic state) with respect to all the degrees 
of freedom but the height of the apical carbon atom, 
and thus different geometries for the triplet and for 
the singlet typically result for the same h value. The 
differences though are minimal as the height of the 
apical carbon atom is the main geometrical parame- 
ter controlling spin alignment in this system, hence 
graphs such as those of Fig. [3] are also representative 
of vertical energy differences. For instance, the pen- 
tagon CC bond length is 2.035 A in the triplet equi- 
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Figure 4: Variation of the open-shell character of the sin- 
glet state as the apical carbon atom moves out of plane. 
The symbols give the weight (c\ + c^) x 100 of the coef- 
ficients in the minimal CAS (2,2) wavefunction with dia- 
batized orbitals described in the main text. The insets 
show isosurfaces at \</>\ = 0.015 A _3//2 of the correspond- 
ing cr-like CASSCF orbitals for representative values of 
h = 0.0, 0.32, 0.68 A, from left to right, respectively. 

librium configuration and increases to 2.081 A in the 
singlet minima, to be compared with dec = 2.007 A 
for the magnetization-free planar structure (Table [i]) 
and dec — 2.467 A in pristine graphene. 

B. Wavefunction calculations 

To analyze in detail the electronic structure in the 
neighborhoods of the vacancy we performed addi- 
tional wavefunction calculations, as described in Sec- 
tion [TTTJ on the geometries obtained at the DFT level, 
fit to the cluster model of Fig. [2] The use of a cluster 
introduces a number of issues related to the finiteness 
of the model (e.g. dependence on the system size, 
discretization of the continuum of states, etc. ) but 
removes some oddities of the supercell approach, e.g. 
the periodic repetition of defects in the same sublat- 
tice which would favor ferromagnetic alignment. We 
first checked that the singlet had the correct open- 
shell character at the planar geometry, namely that 
the wavefunction reads approximately as 

*f =0 oc \...<%<fi\ ± l-^CI 

where <p a and (j)^ are a-like and 7r-like orbitals on the 
apical carbon atoms, respectively, |..| is a shorthand 
for a Slater determinant and the plus (minus) sign 
holds for the triplet (singlet) state. For non-planar 
geometries a— like and tt— like orbitals get generally 
mixed, and the singlet displays both open- and closed- 
shell character. Only if the first dominates the singlet 
can be considered to be the "same" electronic state of 
the triplet but with one spin flipped, and the singlet- 
triplet energy separation is meaningful of an exchange 
coupling. 




Figure 5: CASPT2 energies as function of the displace- 
ment (h) of the apical carbon atom out of the surface 
plane, for the singlet (open symbols) and the triplet (filled 
symbols) states. Also reported for comparison the magne- 
tization constrained DFT results of Fig(3] (dashed lines). 
Energies are referenced to the triplet minimum at the cor- 
responding theory level. The inset shows a blow-up of the 
h ~ region for the singlet. 

To check this, we exploited the invariance of the 
CASSCF wavefunction with respect to rotations of 
the active orbitals, and choose orbitals which maxi- 
mize the overlap (while keeping orthogonality) with 
the above (j) a and <j)^ states of the planar case. In this 
case 

*f=° « Cl |...C^I+C 2 |...^^|+C 3 |...C^I+C4|...C^I 

where c\ and C2 = —c\ represent the open-shell con- 
tribution, and cs and C4 account for the closed-shell 
character (c% = C4 = and C2 = c\ in the triplet). 
In Fig. [4] we report the evolution of the weight 
l c i| 2 + I C2 1 in the normalized wavefunction, as a mea- 
sure of the open-shell character in singlet state. They 
were obtained from simple CAS (2,2) calculations on 
the singlet-optimized geometries -i.e. using just the 
four-determinant wavefunction described above- but 
similar results are found for more elaborate functions. 
Evidently, the system is of open-shell type for a wide 
range of h values, comprising the equilibrium one. 
Only for very large values of h the system prefers a 
closed-shell configuration with "magnetic" properties 
turned off; correspondingly the triplet is pushed much 
higher in energy. 

Multi-configuration SCF wavefunctions obtained 
distributing 4-electrons in 4-orbitals (CAS (4,4)) were 
then optimized for several geometries sampled from 
the singlet and the triplet curves in Fig. [3j and used 
as references for perturbative (CASPT2) calculations. 
The results are shown in Fig. [5] , together with DFT 
ones for comparison, for several values of the h coor- 
dinate, referenced to the triplet minimum. 

It is evident from Fig. [5] that CASPT2 and DFT 
results closely parallel each other for the triplet but 
differ substantially in the singlet. In the latter case, 
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a cusp (due to a likely interaction with higher lying 
electronic states) is only present at the DFT level, 
and smooth out at the CASPT2 level, thereby sig- 
nalling the presence of an avoided crossing. In fact, 
the CAS (4,4) space is sufficiently large to allow us to 
properly describe a number of quasi-degenerate states, 
i.e. those obtained by placing all the four unpaired 
electrons of the vacancy in low-lying states. 

As expected, ferromagnetic coupling is preferred for 
most values of the height of the carbon above the sur- 
face, and a crossing results at about h = 0.5 A; for 
larger values of h, the gain in hybridization energy 
overcomes Coulomb repulsion, and the system show 
increased close-shell character. A minimum occurs in 
the singlet at h = h = 0.38 A (h « 0.4 A at the DFT 
level), and is actually a double minimum, h = =L/i , on 
account of the meaning of the h coordinate. The sin- 
glet is thus a symmetric bistable system with a barrier 
height of E b = 0.09 eV. 

Exchange (Hund) coupling was obtained by the ver- 
tical singlet-triplet energy separation, Jh = AE$T: 
using the geometries optimized for the triplet, and is 
shown in Figj6]as a function of the angle subtended 
by the a dangling bond and the graphene plane. The 
results closely parallel those reported in Fig. [5] since, 
as observed above, the difference between singlet and 
triplet geometries are minimal. The coupling ranges 
from Jh ~ 0.25 eV in the planar configuration to 
about Jh ~ 0.1 eV in the equilibrium configuration 
of the singlet, and thus Jh ~ 0.25 — 0.20 eV seems to 
be appropriate for the ground-state system. 

It is worth stressing at this point, that Jh defined in 
this way is the Hund coupling constant related to the 
geometry-dependent a-like and 7r-like orbitals host- 
ing the unpaired electrons. Its value in the planar 
structure, J H = Jh(6 = 0), gives the Hund cou- 
pling constant in the Anderson impurity model for 
the vacancy^, while its dependence on 6 (at small 
angles) simply reflects the behavior of the hybridiza- 
tion strength^ V CT7r = y/2tg(6) - 2tg 2 (6)) /3 Ae sp 
(Ae sp being the carbon s — p splitting), as confirmed 
by the dashed line in FigjHJ 

Furthermore, despite the limited size (and discrete- 
ness of the energy spectrum) we found no indication 
that the 7r-midgap state is only marginally occupied, 
thereby suggesting that the limit U a7T <C J#/4 applies 
in the above mentioned Anderson model (U a7r is the 
Coloumb repulsion between electrons in the a and in 
the tt midgap state). However, there exists the pos- 
sibility that increasing the cluster size the situation 
reverts and the 7r-midgap state depopulates in favour 
of some other low- lying 7r-state (J = 1/2 problem). 



V. DISCUSSION 

Computed exchange coupling constants are clearly 
too large to have a decoupled response from the two 
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Figure 6: Variation of the Hund coupling constant as a 
function of the angle subtended by the a dangling bond 
and the graphene plane (filled symbols). Dashed line 
is a low-angle weighted-fit of the data to Jh = Jh ~ 
4Atg 2 (0) (1 - 2tg 2 (0)) /3 , which gives J H = 0.268 eV and 
A = 7.59 eV. 

localized electrons to external magnetic fields. The 
presence of a low-lying singlet at energy A above a 
J-paramagnetic ground-state does affect the magneti- 
zation, but only to the extent it modifies thermal pop- 
ulations, that is introducing a temperature- and field- 
dependent correction factor f(fi,H) = A(Be~^ A + 
A) -1 to the thermally averaged magnetic moment 
(here A = sinh[/3 7 # (j + §)] and B = sinh[^], 
/3 = 1/ksT as usual, H is the magnetic field and 
7 is the relevant gyromagnetic ratio). This factor 
has a distinguishing feature of making the moment no 
longer dependent on the reduced field (3H = H/ksT 
only, but is hardly appreciable for /3A > 1 (i.e. 
T < 2000 K(!) for A - 0.2 eV). Only for ft A <C 1 
this factor transforms the J = 1 ground-state magne- 
tization density into twice that of a J — 1/2 moment. 
In practice, this limit attains only if A is vanishing 
small, since the above (2-electron) "atomic" picture is 
challenged at much lower temperatures b y th ermal ex- 
citations out of/into the tt midgap state 50 . Likewise 
for doping which can be used to tune the tt state pop- 
ulation, as it has been recently shown by molecular 
adsorption^. 

All this suggests that the bare vacancy in free- 
standing graphene at low temperature should display 
a J = 1 paramagnetic response and results reported 
by Ref. 8 are consistent with this picture. Decoupled 
a and tt moments, as those observed under better- 
controlled conditions by Nair et aZP^, though, are 
still plausible since the apical carbon atom may be 
easily stabilized out of the surface plane (at about 
h ~ 0.3 A where Jh ~ 0), for instance in the presence 
of a (weakly-binding) substrate^! or ripples. 

With the same token, spin-half residual moments 
may arise because of interaction with foreign species. 
Vacancies are highly reactive species which easily sat- 
urate their dangling a bond and leave a tt magnetic 
moment only. We have checked this by considering 
adsorption of a single hydrogen atom, and found that 
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such a process is both thermodynamically and kineti- 
cally favored at any temperature: no barrier is indeed 
found in the adsorption profile and the binding en- 
ergy is about twice that found for adsorption on a 
bulk site. A deeper analysis of single and multiple hy- 
drogenation comprising static and dynamical aspects 
will be presented shortly. 

Screening of the magnetic impurity by 7r— band 
states (Kondo effect) is a more subtle issue. DFT 
is not able to handle such highly correlated situa- 
tions, and the finite model adopted for the wavefunc- 
tion calculations, along with the limited excitations 
included in the wavefunction, prevent observation of 
any pairing between the "impurity" and the tt band 
states. In the finite size model, such pairing would be 
signalled by the presence of singly excited configura- 
tions where tt states singlet-couple with the impurity 
7r-midgap or a state. Our wavefunction does include 
a number of excitations out of the occupied states in 
the CAS space and, perturbatively, excitations of the 
"core" states at the CASPT2 level, but a detailed anal- 
ysis of the wavefunction such that presented above at 
the CASSCF level is out of question. In the overall 
triplet state, Kondo singlet-pairing would be signalled 
by an increasing derealization of the spin-density 
when enlarging the cluster size but computational 
cost becomes prohibitive to check this [for the sin- 
glet the spin-density vanishes in any case on account 
of the vectorial character of this quantity]. Notice 
though that dynamical mean-field theory with local 
interactions showed no evidence of quenching of the 
7T— related local magnetic moment 18 , in accordance 
with observations of spin-^ paramagnetism in fluo- 
rinated graphene^. In the case of vacancies, screen- 
ing of the <j moment is expected only for non planar 
geometries and, if any, is not compatible with Curie- 
law behavior observed in Ref.s \7\27\ Metallic Kondo 
screening of spin-| impurities has been used to ex- 
plain transport measurements in irradiated graphene 
at different doping levels^, including charge neutral- 
ity, but the interpretation has been questione d 46 ! 47 * , 



and the observed logarithmic increase of resistivity at 
low temperatures related instead to electron-electron 
interactions in the disordered systerrl 4 ^ 4 ^. 



VI. CONCLUSIONS 

We reported on a detailed analysis of the electronic 
and geometric changes that occur upon vacancy for- 
mation in graphene, using both DFT and a high-level 
quantum chemistry method (CASPT2) to overcome 
known limitations of DFT. The picture that emerges 
is that two local magnetic moments coupled to each 
other to give a triplet ground- state, in accordance 
with a report of sp in-1 paramagnetic specie^. Spin- 
half paramagnetism^^, though, can arise in many in- 
stances. Vacancies are highly reactive and easily sat- 
urate their a dangling bond in the presence of foreign 
species. Also, ripples or (weak) coupling to a substrate 
may stabilize a non-planar configuration of the api- 
cal carbon atom, thereby reducing the effective Hund 
coupling constant of the two-electron system and de- 
couple the corresponding local moments. This is the 
likely source of spin-half paramagnetic behavior ob- 
served in Ref.s 17)27] where doping has been shown to 
effectively halve the magnetization density. 

We could not deal with the possible pairing of 
the magnetic moment with the "conduction" 7r-band 
states, because of the limitation of DFT on one hand 
and the use of a finite cluster model (and limited ex- 
citation in the wavefunction) on the other hand. At 
the above level of theory we do not have indication, 
however, of such a coupling. This is consistent with 
the absence of anomalies in the measured susceptibil- 
ity of Ref.s I7)27[ also at finite densities, and in that 
computed (for 7r-moments only) with dynamical mean 
field theory in the presence of local interactions^, and 
suggests that further investigation on the transport 
data measured by Ref. [10] is required for the Kondo 
effect in graphene to be unambiguously identified. 
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